function [inspDist] = calcinspiraldist(m1,m2,f,Sh,f0,fCut);
%
% CALCINSPIRALDIST - calculate the distance to an optimally-oriented binary
% which produces an optimal SNR of 1. 
%
% usage: [inspDist] = calcinspiraldist(m1,m2,f,Sh,f0)
%        [inspDist] = calcinspiraldist(m1,m2,f,Sh,f0,fCut)
% 
% m1,m2    : component masses (in units of M_solar)
% f        ; vector of frequencies at which the PSD is estimated
% Sh       : one-sided PSD of the data (in units of h^2/Hz)
% fUpper   : upper frequency cuf-off. If this is not specified, f_lso is
%            used. 
% inspDist : effective distance to the binary (in MPc)
% 
% P. Ajith, 27.06.2005
% 
% $Id: calcinspiraldist.m 85 2010-01-26 02:01:10Z anand.sengupta $


%% set constants
msolartime  = 4.92579497077314e-06;    % 1 solar mass in geometrical units
timeMparsec = 1.0292712503e14;  % 1 MegaParsec in geometrical units
snr         = 1;                % optimal snr
v_lso       = 1./sqrt(6);       % velocity at the last-stable circular orbit
m           = (m1+m2)*msolartime;
eta         = m1*m2/power(m1+m2,2);
f_lso       = power(v_lso,3)/(pi*m);

if nargin < 6
    fCut = f_lso;
end

delta_f     = f(2)-f(1);                       % freq ewsolution of the psd
N_0         = ceil(f0/delta_f);   % freq bin corresponding to the low-freq cut-off.

%% calculate the effective distance to the prototype binaries. 
N_lso    = ceil(fCut/delta_f);
F_pow    = power(f(N_0:N_lso),-7/3);
S_inv    = 1./Sh(N_0:N_lso);
X        = dot(F_pow,S_inv);
inspDist = (power(m,5/6)/(snr*power(pi,2/3))) * sqrt(5*eta/6) * sqrt(delta_f*X);

%% convert from geometrised units to physical units.
inspDist = inspDist/timeMparsec;



